Many-Body meets QM/MM: Application to indole in water solution 
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Spectral properties of chromophores are used to probe complex biological processes in vitro and in 
vivo, yet how the environment tunes their optical properties is far from being fully understood. Here 
we present a method to calculate such properties on large scale systems, like biologically relevant 
molecules in aqueous solution. Our approach is based on many body perturbation theory combined 
with quantum-mechanics/molecular-mechanics (QM/MM) approach. We show here how to include 
quasi-particle and excitonic effects for the calculation of optical absorption spectra in a QM/MM 
scheme. We apply this scheme, together with the well established TDDFT approach, to indole in 
water solution. Our calculations show that the solvent induces a redshift in the main spectral peak 
of indole, in quantitative agreement with the experiments and point to the importance of performing 
averages over molecular dynamics configurations for calculating optical properties. 

PACS numbers: 



Optical properties of aromatic chromophores embody 
a key facet of cell biology, allowing for a precise interro- 
gation of a variety of biochemical events, including sig- 
naling, metabolism and aberrant processes. These range 
from probing transient interactions between biomolecules 
(proteins and nucleic acids), to protein dynamics and 
fibrillation and plaque formation in neurodegenerative 
diseases. Understanding how the environment tunes 
such optical properties is therefore crucial in structural 
genomics, yet this information is so far mostly lack- 
ing. A powerful tool to address this issue is given by 
the so-called quantum- mechanics/molecular-mechanics 
(QM/MM) methods [HiiilliQ . In this approach, the 
aromatic moiety is treated at quantum mechanical level, 
whilst the environment is described with an effective po- 
tential: the influence of the MM (presumably very com- 
plex and very large) environment is basically included as 
an external potential, and, in case the molecule is cova- 
lently bound to MM region, by a mechanical coupling 
with the environment. 

Most often the QM approach is solved within Density 
Functional Theory (DFT) HQ to study ground state 
properties, and time-dependent DFT (TDDFT) 
when excited states are involved as in the case of the op- 
tical properties TDDFT is computationally 
very efficient, yet its predictive power depends dramati- 
cally on the system and on the functional used to repro- 
duce the exchange and correlation interactions. 

Several approaches, including Post-Hartree-Fock ones 
(configuration interaction and similar methods), 
have been already used to predict optical properties of 
biomolecules. Quantum-many-body techniques (MBPT) 
p^firj . are an attractive alternative, although of course 
they come with a much higher computational cost than 
TDDFT. Strikingly, however, biophysical applications of 
one of the most widely used scheme, the combination of 
the GW method [3 with the Bethe-Salpeter Equation 
(BSE) [13] are so far lacking. The GW method is used 
for the evaluation of the single quasiparticle energies, and 



the BSE to introduce excitonic effects. Keeping in mind 
future biological applications, it is imperative to assess 
the accuracy of a MBPT/MM approach versus the more 
conventional TDDFT/MM one. 

The main assumption in interfacing a QM/MM 
method with TDDFT or MBPT approaches is that the 
optical properties of the chromophore do not involve the 
MM part's electronic structure. Hence, special care has 
to be devoted to the choice of the two regions. 

Here we present MBPT/MM calculations on the in- 
dole ring of the Tryptophan protein residue (Fig. [T]). 
This system appears ideal for such an approach in sev- 
eral respects. First, it is very relevant biologically, as the 
indole ring has been exploited as a spectroscopic tool to 
monitor changes in proteins and to yield information 
about local structure and dynamics. In fact, its spectral 
signatures allow it to be used as a structural probe in 
proteins. Second, it contains a relatively small number 
of atoms (16), which can still be treated at the GW- 
BSE level. Next, the optical gap of liquid water (7 eV 
is larger than the gap of the indole molecule 
(4.3 eV |22|). Under 7 eV the spectra of indole and water 
do not overlap, and it is justified to treat the solvent in a 
classical scheme. Finally, CASPT2 calculations [l^] and 
experimental data [2^ are available, and allow to com- 
pare the changes of the optical properties upon passing 
from the gas phase to aqueous solution. 

We performed QM/MM Car-Parrinello [H simu- 
lations of indole in water by the fully Hamiltonian 
QM/MM scheme 0. Such scheme has been applied to 
a variety of biological systems [2^ . The biomolecule was 
treated, at this step, at the DFT level whilst the solvent 
with the Amber force field [2^. The approach allows for 
an explicit treatment of solvation, in contrast to previous 
studies millll^. 

Indole single quasiparticle energies have been then 
evaluated at the GW level for several snapshots. Fi- 
nally, we solved the BSE to calculate the average absorp- 
tion spectrum and compared the results with the ones 
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obtained within TDDFT. We calculated the indole ab- 
sorbance in water as well as in gas phase. The shift in 
the spectra gives the solvatochromism. 

The GW approximation consists in setting the vertex 
in the Hedin equations [2^ equal to a delta function. 
Under this condition, the time-Fourier transform of the 
proper exchange-correlation self energy, E(r, r',w), is a 
convolution of the Green function G(r, r',a;), with the 
screened Coulomb potential W{r,r' ,uj). The electronic 
bands are obtained by solving the following eigenprob- 
lem: 
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This expression is derived from the Dyson equation in 
Lehmann representation . Vh is the Hartree potential 
of the QM part, [/Q*^ is the electron-ion potential of the 
QM part, while [/QWmm ^-^^ potential felt by the elec- 
trons due to the point charges of the MM part. Finally, 
ej^ are the quasi-particle eigenvalues. Eq. H]) has the 
same form as the KS equation in the presence of an 
external electric field, where the exchange-correlation po- 
tential T4c(r) is replaced by the self energy E(r,r',e^^) 
which acts as a non-local, energy-dependent potential. 
Therefore, the eigenvalue problem described above can 
be solved perturbatively considering the KS equation as 
an unperturbed Hamiltonian and S — Vxc as a perturba- 
tivc term. The quasi-particle eigenvalues are obtained in 
first order approximation: 
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All the Coulomb interactions, and hence also the one 
induced by the classical region, are included in the 
KS eigenvalues ef^, and eigenvectors \4>f^)- In this 
GW/MM scheme we neglect the contribution of the clas- 
sical atoms to the S operator in the same way as it is 
neglected for Vxc in the DFT/MM scheme. 

As a bonus from GW calculations we obtain also the 
four-point independent quasi-particle polarizability re- 
quired to solve the BSE and which reads, in transition 
space jlTj : 
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where /„ is the occupation numbers of level n: the 
quantum plus classical external potential is not explic- 
itly present in the BSE, but indirectly determines all the 
ingredients. 

We performed a 20 ps hybrid QM/MM Car- Parr inello 
simulation of an indole molecule (QM part) [2^ sur- 
rounded by 2000 water molecules treated classically (MM 



FIG. 1: (color online) Indole in water solution. Colors corre- 
spond to the following atomic species: BLUE=N, CYAN=C, 
WHITE=H, RED=0. 



part), with the Amber force field [37|. Such a large num- 
ber of water molecules was necessary to correctly repro- 
duce the physical properties of a disordered system such 
as liquid water at room temperature (300 K). We also 
verified that such a large number of water molecules used 
for the dynamics was also necessary to calculate the op- 
tical absorption spectrum. This became clear already 
at the DFT- independent particle approach (DFT-IPA), 
as shown in Fig[21 we calculated the absorption spec- 
trum for indole with no water (with indole in the same 
"distorted" geometrical configuration as if in water), for 
indole with 2 water molecules, and for indole with 2000 
water molecules. The three spectra are all different. This 
is a strong indication that also for the absorption spec- 
trum it is necessary to include many water molecules in 
the theoretical simulation. This is due to the long range 
electrostatic potential of water which acts on indole. 

Next, we tested our assumption that the solvent can be 
treated classically in the calculation of absorption spec- 
tra by performing TDDFT calculations [5^ for a sys- 
tem where two water molecules were treated at quantum 
level, and the remaining 1998 classically. The position of 
absorption peak is the same (within 0.02 eV) as in the 
case where all waters were treated within MM. Our re- 
sult supports the use of this approach for solutes in water 



For ten snapshots of the QM/MM dynamics (one ev- 
ery two ps) we computed the optical spectra at the in- 
dependent particle level (DFT-IPA) and within TDDFT. 
The final spectrum was obtained by an average over the 
snapshots. The convergence of the spectrum was reached 
already after 6 snapshots, hence the subsequent GW and 
BSE calculations have been performed on only 6 snap- 
shots. 



The calculated DFT and GW HOMO-LUMO gaps 
[ist . averaged over the QM/MM configurations, are 3.8 
eV (with standard deviation a = ± 0.1 cV) and 7.2(cr ~ 
± 0.2) eV, respectively. The GW corrections to the gap 
turned out to be practically constant in all the snapshots 
considered (3.4 ± 0.1 eV). This fact, already found for 
liquid water [2l| . confirms that one can strongly reduce 
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FIG. 2: (color online) DFT-IPA optical spectra of indole in 
2000 water molecules (blue solid line), in 2 water molecules 
(red circles) and without water (black dashed line). 
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FIG. 3: (color online) BSE (a) and TDDFT (b) spectra of 
indole in water. The tiny red dashed lines are the spectrum 
of each snapshot. The red solid line is obtained by an average 
over these spectra. The black dashed line is for indole in vapor 
phase. 

We finally calculate the low energy range of the opti- 
cal spectrum of indole, by GW-BSE and TDDFT, always 
as result of an average over the QM/MM snapshots. In 
FiglHlwe report our results together with the calculated 
absorption spectrum in gas phase. We notice that, in 
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FIG. 4: (color online) BSE optical spectra. Solid blue line: 
Indole in vapor phase. Black dashed line: indole without wa- 
ter molecules, with distorted geometry taken from a snapshot 
corresponding to 13.08 ps of the dynamics. Circles: indole in 
water, spectrum calculated for the same snapshot. 



both approaches the most intense peak (^La in the ex- 
periment) is red-shifted on passing from gas phase to the 
water solution. This agrees with experiments Same 
conclusion was obtained by previous theoretical calcula- 
tions of indole in water based on CASSCF method and 
CASPT2 (2^. In these approaches, the solvent was simu- 
lated by a continuum model with a cavity containing the 
indole molecule. The value we calculate for such a red- 
shift is ~ 0.2 eV, in fair agreement with experiment (0.18 
eV) m. On the contrary the CASSCF and CASPT2 
prediction for the solvent shift is about 0.06 eV only. 
Such an underestimation may depend on the geometri- 
cal distortion of indole molecule caused by temperature 
effects due to the solvent and by an explicit H-bonding 
between water molecules, which were not considered ex- 
plicitly therein. To quantify the effect of the geometry 
distortion on such shift, we performed calculations of in- 
dole switching on and off the water field in order to sepa- 
rate the geometry effect from the electrostatic ones. The 
results are presented for a single snapshot in GW-I-BSE 
(Fig. |4]) and are obtained by performing additional calcu- 
lations for the same QM/MM configuration without the 
water field. The corresponding solvent-shift goes from 
-0.1 eV with water field to -1-0.2 eV (hence, a blue-shift) 
without water field. This emphasises the importance of 
taking into account explicitly the electrostatic interac- 
tion with the solvent, since the geometry distortion alone 
would give, at least for this snapshot, a wrong sign. 

In addition, TDDFT underestimates the energy of the 
^La peak both in gas phase and in solution by ~0.4 eV, 
and BSE-GW overestimates them by ~0.3 cV. As ex- 
pected [i^l, CASSCF is much worse, it overestimates by 
~1 eV or more, whilst CASPT2 is more accurate (^0.13 
cV or less). 

In conclusion, we have included many-body pcrturba- 
tive techniques in a QM/MM scheme. Wc have applied 
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it, together with a TDDFT/MM approach, to study the 
optical properties of indole in water solution. Both meth- 
ods reproduce quantitatively the redshift induced by the 
solvent. Hence, the GW-BSE method could be applied 
to biomolecules in aqueous solution (i.e. in laboratory- 
realizable conditions) in cases where the TDLDA/GGA 
approach docs not work [45ll46j . Our GW-BSE calcu- 
lations further show that the solvent shift is a conse- 
quence of the combination of two effects: the geomet- 
rical distortion of indole molecule in the solvent and the 
electrostatic interaction with the water molecules elec- 
tric dipoles. Both effects, and their sum, depend on the 



particular configuration of the system; this emphasizes 
the need of more than one snapshot (several, indeed) for 
carrying out accurate optical calculations. 

This work opens the way to further applications in 
other bio-relevant molecules, such as proteins and cell 
membranes, for which the evaluation of the optical shift 
enables to understand the nature of their environment. 
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